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A strongly coupled plasma of quark and gluon quasiparticles at temperatures from 1.1T C to 3T C is studied 
by path integral Monte Carlo simulations. This method extends previous classical nonrelativistic simulations 
based on a color Coulomb interaction to the quantum regime. We present the equation of state and find good 
agreement with lattice results. Further, pair distribution functions and color correlation functions are computed 
indicating strong correlations and liquid-like behavior. 
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1 Introduction 

Determining the properties of deconfined quark-gluon plasma (QGP) is one of the main challenges of strong- 
interaction physics both theoretical and experimental. Many features of this matter were experimentally dis- 
covered at the Relativistic Heavy Ion Collider (RHIC) at Brookhaven. The most striking result, obtained from 
analysis of these experimental data [1], is that the deconfined quark-gluon matter behaves as almost perfect fluid 
rather than as a perfect gas, as it could be expected from the asymptotic freedom. 

From the theory side, the most fundamental way to compute properties of strongly interacting matter is pro- 
vided by lattice QCD, see the recent review [2]. Interpretation of these computations requires application of 
various QCD motivated, albeit schematic, models simulating various aspects of the full theory and allowing for 
a deeper physical understanding. The above mentioned strongly correlated behavior of the QGP is expected to 
show up in long-ranged spatial correlations of quarks and gluons which, in fact, may give rise to liquid-like and, 
possibly, solid-like structures. This expectation is based on very similar behavior observed in electrodynamic 
plasmas, as was pointed out e.g. in Refs. [3,4]. This similarity has been exploited to formulate classical non- 
relativistic models of a color Coulomb interacting QGP [4] which is numerically analyzed by classical molecular 
dynamics simulations. While this has allowed to incorporate nonideality effects, quantum effects were either 
neglected or included phenomenologically via a short-range repulsive correction to the pair potential, e.g. [4]. 
However, such a rough model may become a critical issue at higher densities where quantum and Fermi statistics 
effects of the quarks should have a strong influence on the properties of the QGP. Similar models had been used 
in electrodynamic plasmas and showed poor behavior in the region of strong wave function overlap, in particular 
at the Mott density. 

This difficulty can be eliminated by deriving effective quantum potentials, as was shown by some of the present 
authors before [5-7]. Following an idea of Kelbg [8] quantum corrections to the pair potential can be rigorously 
derived in perturbation theory with respect to the coupling parameter [9]. To extend the method of quantum 
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potentials to larger coupling an "improved Kelbg potential" was derived which contains a single free parameter 
which is fitted to the exact solution of the quantum-mechanical two-body problem and exactly reproduces the 
thermodynamic properties up to moderate couplings [6]. However, this approach may fail if the system forms 
bound states of more than two particles leading to a break-down of the pair approximation for the density matrix, 
as was demonstrated in Ref. [6]. A superior approach which does not have these limitations, is to use the original 
Kelbg potential in path integral Monte Carlo (PIMC) simulations which effectively map the problem onto a high- 
temperature weakly coupled and weakly degenerate one. This allows one to rigorously extend the analysis to 
strong couplings and is, therefore, the method of choice for the present purpose. Beside the nonideality and 
quantum effects our approach takes into account the effects of the Fermi (Bose) statistics of quarks (gluons) by 
a proper antisymmetrization (symmetrization) of the N— body density matrix. For temperature and density of 
the QGP considered in this paper which are similar to the conditions of Ref. [4] these effects are very important 
as the quasiparticle thermal wave length is larger than the average interparticle distance. This is, in particular, 
important for the behavior of the pair distribution functions (see below). 

Here, we develop a PIMC approach to the strongly coupled QGP which takes the Fermi (Bose) statistics of 
quarks (gluons) and quantum degeneracy selfconsistently into account. This method has been successfully ap- 
plied to strongly coupled electrodynamic plasmas before, e.g. [10-12]. Examples are partially ionized dense 
hydrogen plasmas where liquid-like and crystalline behavior was observed [13, 14]. Moreover, also partial ion- 
ization effects and pressure ionization could be studied from first principles [15]. The same methods have been 
applied also to electron-hole plasmas in semiconductors [16, 17], including excitonic bound states, which have 
many similarities to the QGP due to the smaller mass differences as compared to electron-ion plasmas. 

In this paper we present first exploratory PIMC simulations of a nonideal quark-gluon plasma. The main goal 
is to test this approach for ability to reproduce the equation of state known from lattice data [2]. To this end we 
use the simplest model of a QGP consisting of quarks, antiquarks and gluons interacting via a color Coulomb 
potential due to Gelman et al. [4] with several approximations for the temperature dependence of the quasiparticle 
masses. We report surprisingly good agreement with the lattice data for one of the parameter sets, which gives us 
confidence that the model correctly captures main properties of the nonideal QGP. 

The paper is organized as follows. In Sec. [2] we introduce the model and approximations used which is 
followed by an overview on our PIMC simulations in Sec. [3] Sec.|4]contains our results on the equation of state 
and on the various pair distribution functions of the QGP, and we conclude in Sec. [5] with a discussion of the 
results together with an outlook on further improvements of the approach. 

2 Theoretical Model 

Our model is based on precisely the same assumptions as those in Ref. [4] which are summarized as follows: 

I: All particles (quarks and gluons) are heavy, i.e., m > T, where m is the mass of a particle and T the 
temperature and, therefore, they move non-relativistically. This assumption is based on the analysis of 
lattice data [18,19]. 

II: Since the order of magnitude of quark and gluon masses, deduced from the lattice data [18, 19] is the same, 
we do not distinguish these masses and put them equal. Moreover, because of the latter we do not distinguish 
between quark flavors. 

Ill: The interparticle interaction is dominated by a color-electric Coulomb interaction, see Eq. ([TJ. Magnetic 
effects are neglected as sub-leading ones, in the nonrelavistic limit. 

IV: The color operators t a are substituted by their average values, i.e. by classical color vectors, relying on the 
fact that the color representations are large. 

The quality of these approximations and their limitations were discussed in Ref. [4]. 

Aiming at a first test of this model in PIMC simulations, in this paper, we are going to consider the QGP only 
at zero baryon density. Therefore, this model requires the following quantities as an input: 

1. the temperature dependence of the quasiparticle mass, m(T), 
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2. the density of particles, n(T), at a given temperature, following [4], where we assume that the numbers of 
quarks, antiquaries and gluons are approximately equal, 

3. the coupling constant, g 2 (T), at a given temperature, see Eq. ([TJ. [Note that, because of the running coupling 
in the QCD, g 2 generally depends on T]. 

All the input quantities should be deduced from the lattice data or from an appropriate model simulating these 
data. 



3 Path integral Monte Carlo Simulations 

As discussed in section [2] we consider a three-component QGP consisting of N = N q + N q + N g quasipar- 
ticles, where N q ,N q ,N g are, respectively, the number of (dressed) quarks, antiquarks and gluons in thermal 
equilibrium, so the temperature-dependent Hamiltonian can be written as Hp = Kp + Ug, where ft = 1/ksT, 
is the inverse temperature and ks is Boltzmann's constant. Here we introduced the kinetic and color Coulomb 
interaction energy of the qu as ip articles 

ft -ST Pt fro _ 1 V^X^ Cpt9 2 {Q P \Qt) , n 



Here the Q p denote Wong's color variables which are 3D unit vectors, and the constants C p t = Ct p are products 
of eigenvalues of the Casimir operator [4]: C qq — C qq = C qq = 4/3, C qg — C qg = 2 and C gg = 3 with 
.9 2 = 27r. ' ' 

The thermodynamic properties in the canonical ensemble with given temperature T and fixed volume V are 
fully described by the density operator p = e~@ H jZ with the partition function (normalization constant) 

Z(N q , N q , N g , V; p) = - ^ m j J drdQp(r, Q, a; 0), (2) 

1 1 9 a y 

where p(r, Q, a; /3) denotes the diagonal matrix elements of the density operator at a given value a of the total 
spin. In Eq. ([2]), r — {r q ,r q ,r g } and Q = {Q q ,Q q ,Q g } are the spatial and color coordinates, while a = 
{cr q , a q } are the spin degrees of freedom, i.e. r a = {r lya . . . n, a • ■ ■ r Na ,a}, Qa = {Qi,a ■ ■ ■ Ql,a ■ ■ ■ QN a ,a} and 
<y a = {&i,a ■ • • &la ■ ■ ■ CN a ,a} with a= q,q, g. 

In order to calculate thermodynamic functions, the logarithm of the partition function has to be differentiated 
with respect to thermodynamic variables. For example, for pressure and internal energy follows 

/3P = dinZ/dV = [a/3Vd\nZ/da] a=1 , (3) 
0E = -pdhxZ/dp, (4) 

where a = L/L is a length scaling parameter. 

Of course, the exact density matrix of interacting quantum systems is not known (particularly for low temper- 
atures and high densities), but it can be constructed using a path integral approach based on the operator identity 
e~^ H — e~ Al3H ■ e~ Al3H . . . e~ A/3H , where the r.h.s. contains n + 1 identical factors with A/3 = /3/(n+ 1), 
which allows us to rewrite the integral in Eq. Q 

J dr<®dQ<® p(q {0) ,Q°, <j; (3) = J dr^dQ^ . . . dr^dQ™ p« • p^ . . . p<»> x 

E E E E(- 1 )^ +Kp - s ^ hhp a °') p q p,p a p {n+ \^ =rm ^ • (5) 

<T P q P, P g 

The spin gives rise to the spin part of the density matrix (S) with exchange effects accounted for by the permu- 
tation operators P q , P q and P g acting on the quark, antiquark and gluon spatial r(" +1 ) and color Q(" +1 ) coordi- 
nates and spin projections a'. The sum is over all permutations with parity np q and /tp r In Eq. (|5| the index k = 
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1 . . . n + 1 labels the off-diagonal high-temperature density matrices p( k > = p ^ r ( fe i)q(* 1 ) ) 7-( fe )Q( fe ); A/3) = 
^ r (fe-i)| e -A/3H| r (fe)^^Q(fc-i) _ Q(k)y Accordingly each particle is represented by a set of n + 1 coordinates 
("beads"), i.e. the whole configuration of the particles is represented by a 3(N q + N q + N g )(n + 1) -dimensional 
vector f = {r^ q , . . . r[ n q +1 \ rf^ q . . . r^ q , . . . r^^jr^ . . . r^^'} and a 2(iV g + A, +A ff ) -dimensional color 

vector Q = {Qf q , Q^g, ■ ■ ■ , g }- The main contributions to the partition function come from configurations 
in which the 'size' of the cloud of beads of quasiparticles is of the order of their thermal wavelength, whereas 
typical distances between beads of each quasiparticle are of the order of the wavelength taken at the (n + 1)- 
times higher temperature A/3. To determine the pressure or total energy in the path integral representation ([5]) 
each high-temperature density matrix has to be differentiated in turn according to expressions Q and Q. 

Let us now consider approximations for the high-temperature density matrices p s k- An approximation which is 
suitable for direct PIMC simulations has the following form, generalizing the electrodynamic plasma results [12] 
to the case of an additional bosonic species (the gluons): 

n N 

Psk(r,P) = C^ 9 e-^-^nn^ dct ll^ ,1 |l-* x P er l^ n,1 |l' (6) 

1=1 p=l 

where s(k) is the number of quarks (antiquarks) with the same spin projection, antisymmetrization and sym- 
metrization are taken into account by the symbols "det" and "per" denoting the determinant and permanent, 
respectively. Further, we introduced the total color interaction energy 

n 

U(r, Q,(3)=Y,U (r« , Q, A/?) /(n + 1), (7) 

Here, the result is rewritten in terms of dimensionless coordinates r which depend on the dimensionless distances 
between neighboring beads, £± , ■ ■ ■ $f , according to r = [r^; + j/ 1 ); A 2 ^ + ?/ 2 ); . . . ; + y^}, 
with y n — AA Y^k=i Further, we introduced the high-temperature De Broglie wavelength, AA^ = 

2" 



2nfi 2 A[3 1 'm a (/3) , and the exchange matrix </>' = exp 



pp 

The path integral representation of the density matrix is exact in the limit n — » oo. For any finite num- 
ber n, the error of the above approximations for the whole product on the r.h.s. of Eq. |5]) is of the or- 
der l/(n + 1) whereas the error of each high-temperature factor is of the order l/(n + l) 2 , as was shown 
in Ref. [12]. Our approximation of the high-temperature density matrix is given by products of two-particle 
density matrices = p$ + 0[(l/n + l) 2 ], where is the kinetic density matrix, while p^ — 
e -Aj3U{r {% 1 ,Q)j^ r (i-i) _ r W^ where U denotes the pair sums of the off-diagonal two-particle effective quan- 
tum potentials. These potentials are straightforward generalizations of the corresponding potentials of electro- 
dynamic plasmas [6, 12] to the case of color Coulomb interaction, <& p *(r p , r^, r t , r' t , Q). In the following we 
will use the diagonal element (r' = r p ,r' t = r t ) which will be called "color Kelbg potential" and depends 
only on a single distance. It is derived by approximating the off -diagonal matrix elements of the effective binary 
interaction by the diagonal ones at the center coordinate Q ab (r, r' ,Qa,Qb'-, A/3) ~ <I> ab ( z: ±^, , Q a , Qi,; A/3) or 
$ ab (r, r'; , Q a , Qb, A/3) « [$ afc (r, , Q a , Q b , A/3) + $ afc (r', Q a , Q b , A/3)]/2. The result for the diagonal color 
Kelbg potential is (we retain the same notation as before) 



(r p , r p , r t , r t , Q p , Q t , A/3) = ( | r pt | , Q p , Q t , A/3) 
C pt g 2 (Q P \Qt) 



4ir\ pt Xpt 



1 - e x * + ^x pt (1 - erf (x pt )) 



(8) 



where x pt = |r p — r t |/A pt . Note that the color Kelbg potential approaches the color Coulomb potential at 
distances larger than the De Broglie wavelength. Most importantly, it is finite at zero distance (it is of the order 
of T), removing in a natural way the classical divergences which makes any artificial cut-offs obsolete. 

Finally let us comment on the treatment of the exchange properties of quarks, antiquarks and gluons. The 
density matrix ([6]) has been transformed to a form which does not contain an explicit sum over permutations and 
thus no sum of terms with alternating sign (in the case of quarks and antiquarks). Instead, the whole exchange 
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problem is contained in exchange matrices from which we have to compute the determinant (for quarks and 
antiquarks) or the permanent (for gluons), 



\\<f>pLi 



J pt,ol\\sk - 

■^\(r P a-n,n)+v^ (Q p \Q t ) + 1) 



«,i (Q P \Qt) + l 
pt 2 



n,l (QM+l 

ol • 2 



a|(ro.i-r, lB )+»f| a (Q \Qi) + l 



(9) 



As a result of the spin summation, the matrix carries subscripts sk denoting the number of quarks and antiquarks 
having the same spin projections. 



4 Numerical results 

In this section we present results of our simulations. Details of our path integral Monte Carlo simulations have 
been discussed before in a variety of papers and review articles, e.g. [21] and references therein, and will not be 
repeated here. The main idea of the simulations consists in constructing a Markov chain of configurations which 
differ by the particle coordinates (including all beads). In addition to the case of electrodynamic plasmas, here we 
also randomly modify the color variable Q of all particles until convergence is achieved. For the results presented 
below we used a cubic simulation box with periodic boundary conditions. The number of particles was equal to 
N = N q + Ng + N g — 40 + 40 + 40 = 120, and the number of high-temperature factors (beads), n = 20. 




Lattice EOS 
PIMC EOS [4] 
PIMC EOS n=3fm 



2.5 T/T 3.0 



Fig. 1 Equation of state of the QGP from PIMC simulations compared to lattice data of Ref. [2]. Line 1 - lattice data [2], 
line 2 - PIMC results for parametrization a, line 3 - PIMC results for parametrization b, see text. 

In Fig. [T] we present results for the QGP equation of state (EOS) obtained from lattice calculations, cf. Line 
1, [2] and from our PIMC calculations according to Eq. ^ with N q = Nq = N g = N/3, based on two different 
parametrizations of the quasiparticle model: 

a: Temperature dependence of quasiparticle density, mass and coupling constant are chosen according to 
m(T)/T c = 0.9/(T/T c - 1) + 3.45 + 0.4T/T C [4, 19]. Results for n{T) = 0.244T 3 and g 2 {T) = 2ir 
are shown by Line 2 in Fig. |T| 

b: Same mass and coupling constant as in a, but with T-independent quasiparticle density, n w 3/m -3 . The 
results are shown by Line 3. 
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As seen from Fig.[T] the version with the T-dependent density (a) results in substantial deviations from the lattice 
data. In contrast, the constant-density version (a) gives a surprisingly good agreement with the lattice results in 
the whole range of temperatures down to values as low as T — l.2T c . Despite the simplicity of model b it seems 
to capture basic trends of the global thermodynamic properties of the QGP. 

With the PIMC simulations we are now able to analyze more in detail additional properties and the internal 
structure of the QGP which can be understood from the pair distribution functions, see below. First we note that 
the QGP in the studied temperature range is, in fact, quantum degenerate. This follows from the degeneracy 
parameter \ a — n a ts? a , where A = h/ \/2irm a kBT is the thermal De Broglie wave length of the quasiparticle 
of species "a" (here it is the same for quarks, and gluons). In the studied temperature interval, \ is practically 
constant and equal 4.1. From this we expect that the finite extension of the quasiparticles is relevant and also spin 
statistics (e.g. the Pauli principle) should play a significant role. At the same time, the relatively moderate value 
of x indicates that the chosen number of high-temperature factors in the PIMC simulations is appropriate. 





Fig. 2 Pair distribution functions (upper panel) and color pair distribution functions (lower panel) of identical particles (left 
column) and different particles (right column) at temperature T/T c = 3 for the parametrization b. The mean interparticle 
distance is (r) ja — 0.06, where the length scale is defined as a — 2ivhc/kBT c = 7.27 fm and T c = 175MeV. 

Let us now consider the spatial arrangement of the quasiparticles in the QGP more in detail. To this end we 
analyze the pair distribution and color pair distribution functions. The pair distribution functions (PDF) g a b(r) 
give the probability to find a pair of particles of type "a" and "b" at a certain distance r. In a non-interacting 
classical system, g a b{r) = 1, whereas interactions and spin effects cause re-distribution of particles. The PDF is 
defined according to 

9ab{Ri,R2) = g ab (Ri-R2) = ^J2[ drd ^^ Rl - r ^ R2 - r ^P^^ a '^^ (10) 

CT V 

Z = Z(N q ,N q ,N g ,V;P)N q \N q \N g \, (11) 
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and results for the PDF at temperature T/T c — 3 are shown in Fig. [4] top panel. Let us first consider the PDF 
between identical particles, see top left figure. At large distances, r/a> 0.15 where a — hc/ksTc — 7.27fm, 
all functions coincide, approaching unity, as in the ideal gas case. However, there is a drastic difference in the 
behavior of the PDF of quarks and gluons (the anti -quark PDF is identical to the quark PDF) and small distances. 
While the gluon PDF increase monotonically when the distance goes to zero, the PDF of quarks (and antiquarks) 
exhibits a broad minimum. This difference can be understood by spin statistics effects. In fact, for the present 
conditions, the De Broglie wavelength A equals approximately O.lcr, i.e. the different behavior of gluons and 
quarks starts to appear at distances smaller than 1.5A. The enhanced population of low distance states of gluons 
is due to bosonic exchange and color Coulomb attraction, see below. In contrast, the depletion of small distance 
states of quarks is a consequence of the Pauli principle. In an ideal Fermi gas g(r) equals zero for particles with 
the same spin projection and one for particles with opposite spin projection, in the limit r — > 0. As a consequence, 
the spin averaged PDF approaches 0.5 [22]. This low-distance behavior is also observed in a nonideal dense 
astrophysical electron-ion plasma and in nonideal a electron-hole plasmas in semiconductors [22,23]. This is 
exactly the value we observe for the quark-quark PDF at small distances which we, therefore, attribute to the 
Pauli principle. An exception is the behavior at very small distances, r < 0.02c. Also, the abrupt increase of 
g around r = O.lcr is a-typical compared to electrodynamic plasmas. We, therefore expect that this behavior is 
caused by the particular properties of the color Coulomb interaction. 

Let us now consider the PDF of different particles, see top right part of Fig. 2. Here all curves show similar 
behavior. At small distances, r < 0.05c, a strong increase is observed which resembles the behavior of the 
gluon-gluon PDF, cf. top left figure. At larger distances, all PDF's equal one. This increase of the PDF at 
small distances is a clear manifestation of an effective pair attraction of quarks and antiquarks as well as quarks 
(antiquarks) and gluons. This may, at first sight, seem surprising because all Casimir indices C p t are positive, see 
above, indicating a repulsive character of the pair potential ([lj. However, this potential still contains the scalar 
product of the color vectors, and the net attraction could be understood if, on average, the color vectors of nearest 
neighbor quasiparticles of any type are anti-parallel. 

This hypothesis is readily verified from the PIMC simulation data. To this end, we define the color pair 
distribution function (CPDF) by generalizing the definition ( fT0] i according to 

c ab (R 1 -R 2 ) = lYl J drdQ{Q a 1 \Ql)8{R 1 -r a 1 )8{R 2 -r\)p{r,Q,a;f3), (12) 
a v 

which is straightforwardly computed during the PIMC simulations together with the traditional PDF. The results 
are shown in the lower panel of Fig. 2. We immediately observe that all CPDF's are negative at small distances 
indicating anti-parallel orientation of the color vectors of all neighboring quarks (antiquarks) and gluons as well as 
quarks and antiquarks, clearly confirming the origin of the effective quasiparticle attraction seen in the functions 
g a b for a ^ b. We now turn to the CPDF of identical particles, see bottom left figure. All functions are non- 
positive everywhere. The minimum of c qq close to r = explains the increase of g qq above the value 0.5 at 
small distances. Most striking is the deep minimum of the gluon CPDF, c gg , at small distances. It again confirms 
the antiparallel arrangement of the color vectors of neighboring gluons whereas the much lower value of the 
minimum, compared to that of the quark CPDF, is due to the absence of the Pauli principle and the larger value 
of the Casimir index C gg compared to C qq . This deep minimum explains the high maximum of the gluon PDF 
9gg- 

Let us summarize the local ordering of the QGP at the temperature T/T c = 3. We observe only weak signs 
of a spatial ordering, cf. the peak of the quark PDF around r = 0.11c, which may be interpreted as emergence 
of liquid-like behavior of the QGP. Much more pronounced is the short range structure of nearest neighbors. 
The QGP lowers its total energy by minimizing the color Coulomb interaction energy via a spontaneous "anti- 
ferromagnetic" ordering of color vectors of gluons. This gives rise to a clustering of gluons which is accompanied 
by a weak tendency of clustering of quark pairs with anti-parallel spins. We also observe clusters of quarks, 
antiquarks and gluons. To verify the relevance of these trends a more refined spin-resolved analysis of the PDF 
and CPDF is necessary, together with simulations in a broader range of temperatures which are presently under 
way. 
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5 Discussion 

Experimental data on the quark-gluon plasma and the hadronization transition give rise to numerous challenges 
to the theory, see, e.g. [1,20] and references therein. Of particular interest is the question why the quark-gluon 
matter behaves as an almost perfect fluid rather than as a perfect gas, as it could be expected from the asymptotic 
freedom. Quantum Monte Carlo simulations based on the quasiparticle picture with color Coulomb interactions 
help us to answer this question. Indeed, the ratio of the potential energy of the system to the kinetic one, obtained 
in these simulations, turns out to be in the range from 1 to 3, depending on the temperature. This certainly 
corresponds to a liquid-like rather than a gas-like behavior. 

We have shown that the PIMC method captures main trends of the equation of state (even near the critical 
temperature) and may also yield valuable insight into the internal structure of the QGP, in particular into the pair 
correlation functions. Our PIMC simulations also allow for a selfconsistent analysis of cluster and bound state 
formation in the QGP. Similar questions have been successfully studied before in dense astrophysical plasmas 
[13] and electron-hole plasmas in semiconductors [15]. In fact, first indications for clustering in the QGP have 
been observed and will be studied in more detail in the future. 

The PIMC method is not able to yield dynamical and transport properties of the QGP. One way to achieve 
this is to develop semiclassical molecular dynamics simulations. In contrast to previous MD simulations where 
quantum effects were included phenomenologically via a short range potential [4] a more systematic approach 
has been developed for electron-ion plasmas [5,6]. There an effective quantum pair potential has been derived 
from quantum Monte Carlo data which should also be possible in application to the QGP. Finally, another very 
promising approach to study the dynamical and transport properties of strongly coupled Coulomb systems is 
based on the Wigner formulation of quantum dynamics [24] which should also be applicable to the quark gluon 
plasma. 

Acknowledgements 

We acknowledge stimulating discussions with Prof. B. Friman and Prof. M.I. Polikarpov and financial support 
by the Innovationsfond Schleswig-Holstein. Y.I. and V.S. were partially supported by the Bundesministerium filr 
Bildung und Forschung (BMBF project RUS 08/038). Y.I. acknowledges support of the Russian Federal Agency 
for Science and Innovations (grant NSh-3004.2008.2). 

References 

[1] E. Shuryak, Prog. Part. Nucl. Phys. 6 2, 48 (2009)., E. Shuryak, Prog. Part. Nucl. Phys. 53, 273 (2004) 

[2] A. Bazavov, et al., |arXiv:0903.4379| /l [hep-lat]. 

[3] M.H. Thoma, IEEE Trans. Plasma Science 32, 738 (2004) 

[4] B.A. Gelman, E.V. Shuryak, and I. Zahed, Phys. Rev. C 74, 044908 (2006); ibid. 74, 044909 (2006) 

[5] A. Filinov, M. Bonitz, and W. Ebeling, J. Phys. A: Math. Gen. 36, 5957 (2003) 

[6] A. Filinov, V. Golubnychiy, M. Bonitz, W. Ebeling, and J.W. Dufty, Phys. Rev. E 70, 04641 1 (2004) 

[7] W. Ebeling, A. Filinov, M. Bonitz, V. Filinov, and T. Pohl J. Phys. A: Math. Gen. 39, 4309 (2006) 

[8] G. Kelbg, Ann. Physik (Leipzig) 12, 219 (1962); 13, 354 (1963) 

[9] The idea to use a Kelbg-type effective potential also for quark matter was proposed independently by K. Dusling, and 
C. Young, arXiv:0707.2068 <i2. However, their potentials are limited to weakly nonideal systems. 

[10] M. Bonitz et al, Phys. Plasmas 15, 055704 (2008) 

[11] M. Bonitz etal. J. Phys. A: Math. Gen. 36, 5921 (2003) 

[12] V.S. Filinov, M. Bonitz, W. Ebeling, and V.E. Fortov, Plasma Phys. Control. Fusion 43, 743 (2001) 

[13] V.S. Filinov, M. Bonitz, and V.E. Fortov, JETP Lett. 72, 245 (2000), [Pisma v ZhETF, 72, 361 (2000)] 

[14] M. Bonitz, V.S. Filinov, V.E. Fortov. P.R. Levashov, and H. Fehske, Phys. Rev. Lett. 95, 235006 (2005) 

[15] V.S. Filinov, M. Bonitz, P.R. Levashov, V.E. Fortov, W. Ebeling, M. Schlanges, and S.W. Koch, J. Phys. A: Math. Gen. 
36, 6069 (2003) 

[16] M. Bonitz, V.S. Filinov, V.E. Fortov. P.R. Levashov, and H. Fehske, J. Phys. A: Math. Gen. 39, 4717 (2006) 

[17] V.S. Filinov, H. Fehske, M. Bonitz, V.E. Fortov, and P.R. Levashov, Phys. Rev. E 75 , 036401 (2007) 

[18] P. Petreczky, F. Karsch, E. Laermann, S. Stickan, and I. Wetzorke, Nucl. Phys. Proc. Suppl. 106, 513 (2002). 

[19] J. Liao and E.V. Shuryak, Phys. Rev. D 73, 014509 (2006) 

[20] T.S. Biro, J. Phys. G: Nucl. Part. Phys. 35, 044056 (2008) 



Copyright line will be provided by the publisher 



cpp header will be provided by the publisher 



9 



[21] A.V. Filinov, and M. Bonitz, in: "Introduction to Computational Methods for Many-Body Physics", M. Bonitz and D. 

Semkat (eds.), Rinton Press, Princeton (2006) 
[22] V.S. Filinov, V.E. Fortov, M. Bonitz, and D. Kremp, Phys. Lett.A 274, 228 (2000) 
[23] V.S. Filinov, H. Fehske, M. Bonitz, V.E. Fortov, P. Levashov, Phys. Rev. E 75, 036401 (2007) 
[24] V.S. Filinov, P. Thomas, I. Vagra, T. Meier, M. Bonitz, V. Fortov, S.W. Koch, Phys. Rev. B. B65, 165124 (2002) 



Copyright line will be provided by the publisher 



